out<-out[[1]]
ideal.points<-grep("theta",colnames(out))
difficulties<-grep("beta",colnames(out))
locations<-grep("alpha",colnames(out))

normalize.model<-function(mat) {
	my.mat<-mat
	ref.mean<-mean(my.mat[ideal.points],na.rm=T)
	ref.sd<-sd(my.mat[ideal.points],na.rm=T)

	my.mat[ideal.points]<- (my.mat[ideal.points] - ref.mean)/ref.sd

	my.mat[locations]<- my.mat[locations] - ref.mean * my.mat[difficulties]

	my.mat[difficulties]<- my.mat[difficulties] * ref.sd

	my.mat
}

my.model<-as.mcmc(t(apply(out,1,normalize.model)))

## Assess convergence
my.judge.geweke<-try(geweke.diag(my.model[,grep("theta",colnames(my.model))]))

## Start getting estimates out
holder<-summary(my.model,quantiles=c(0.025,0.05,0.5,0.95,0.975))
xbar<-holder$statistics[grep("theta",names(holder$statistics[,1])),1]
alphabar<-holder$statistics[grep("alpha",names(holder$statistics[,1])),1]
betabar<-holder$statistics[grep("beta",names(holder$statistics[,1])),1]

## You could compare it to the original matrix
votes<-cast(data,casenum~judge,value="vote")
votes<-votes[,-1]

## Or you could build a small matrix
my.mat<-matrix(nrow=nrow(data),ncol=3)
colnames(my.mat)<-c("alpha","beta","theta")
my.mat[,"alpha"]<-alphabar[case[1:length(data$vote)]]
my.mat[,"beta"]<-betabar[case[1:length(data$vote)]]
my.mat[,"theta"]<-xbar[as.numeric(factor(data$judge))[1:length(data$vote)]]

get.mu<-function(x) {
	mu<- -1 * x[1] + (x[2]*x[3])
	mu
}
mu<-apply(my.mat,1,get.mu)

predprob <- pnorm(mu)
pred <- predprob >= 0.5
correct <- data$vote == pred
correct[which(is.na(data$vote))] <- NA

## Percentage correct
overall.percent <- (sum(correct, 
            na.rm = T)/sum(!is.na(correct))) * 100 ## 
yea.percent <- (sum(correct[data$vote == 
            1], na.rm = T)/sum(!is.na(correct[data$vote == 1]))) * 
            100 ## 
nay.percent <- (sum(correct[data$vote == 0], na.rm = T)/sum(!is.na(correct[data$vote == 
            0]))) * 100 ## 

## Null model
null.pred<-rep(1,length(data$vote))
null.correct <- data$vote==null.pred
null.percent=(sum(null.correct,na.rm=T)/sum(!is.na(null.correct)))*100 ## 

## percentage correct by legislator

lp <- tapply(correct, factor(data$judge), mean) * 100
#names(lp)<-judge.names

##
## logLik
pijt_1<-log(predprob)
pijt_0<-log(1-predprob) ## optionally, pijt_0<-log(1-predprob+1e-7)

pijt_1<-pijt_1*(data$vote==1)
pijt_0<-pijt_0*(data$vote==0)

my.logLik<-sum(pijt_0,pijt_1,na.rm=T)

my.A<-length(which(!is.na(data$vote)))

GMP<-exp(my.logLik/my.A)
GMP

## Now calculate the GMP for the null model
null.predprob<-rep(mean(data$vote),length(predprob))
null.pijt_1<-log(null.predprob)
null.pijt_0<-log(1-null.predprob)

null.pijt_1<-null.pijt_1*(data$vote==1)
null.pijt_0<-null.pijt_0*(data$vote==0)

null.logLik<-sum(null.pijt_0,null.pijt_1,na.rm=T)

null.GMP<-exp(null.logLik/my.A)
null.GMP

## get results
results.df<-data.frame(holder$quantiles[grep("theta",rownames(holder$quantiles)),])
names(results.df)<-c("Lo","X5","Median","X95","Hi")
results.df$Judge<-judge.names
my.order<-order(results.df$Median)
results.df$Judge<-factor(results.df$Judge,levels=results.df$Judge[my.order],ordered=TRUE)

ideal.plot<-ggplot(data=results.df,aes(x=Judge,y=Median,ymax=Hi,ymin=Lo)) + geom_pointrange() + theme_bw() + scale_y_continuous("Position") + scale_x_discrete("Judge") + coord_flip()

